Multifractality: generic property of eigenstates of 2D disordered metals 
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The distribution function of local amplitudes of eigenstates 
of a two-dimensional disordered metal is calculated. Although 
the distribution of comparatively small amplitudes is gov- 
erned by laws similar to those known from the random matrix 
theory, its decay at larger amplitudes is non-universal and 
much slower. This leads to the multifractal behavior of in- 
verse participation numbers at any disorder. From the formal 
point of view, the multifractality originates from non-trivial 
saddle-point solutions of supersymmetric <r-model used in cal- 
culations. 



Fractality of some unusual objects of condensed mat- 
ter physics has been intensively discussed during the last 
decade Q . This notion has also penetrated to the theory 
of disordered systems for characterizing the vicinity of 
the localization transitions. For example, using the coef- 
ficients of the inverse participation numbers (IPN) known 
from the renormalization group treatment in 2 + e dimen- 
sions j| Castellani and Peliti [|] have suggested the mul- 
tifractal structure of wave functions at the localization 
threshold. Later, due to considerable efforts in numeri- 
cal simulations, it has been recognized that, at criticality, 
the wave functions are really multifractal which gives rise 
to a spectrum of critical exponents B . 

The concept of multifractality of wave functions has 
been introduced |3j^|] as a way to characterize their com- 
plexity in the pre-localized regime and originates from 
the numerically observed non-trivial dependencies of the 
inverse participation numbers t n (V) on the volume V of 
a system. By definition, the latter are the moments of the 
distribution function f(t) of local amplitudes t = \tp(r)\ 2 
of wave functions at an arbitrary point r inside a sample, 



tn = Jo°° t n f(t)dt, (1) 



/(*) = iVu)- 1 (J2 S(t - \Mr)\ 2 )S(e - e a )J 

where ip a (r) and e Q are eigenfunctions and eigenvalues 
of state a of a confined system. The symbol () stands 
for averaging over disorder; v is the average density of 
states. The distribution function and the wave functions 
are properly normalized such that to = 1 and t\ = V -1 . 



The inverse participation numbers indicate very sensi- 
tively the degree of disorder- induced-localization of states 
through their dependence t n (V) on the volume of the sys- 
tem. In a pure metal or a ballistic chaotic box where the 
wave functions extend over the whole system one has 



t n oc l/V n . 



(2) 



If disorder makes the localization length L c much shorter 
than the sample size L ~ V 1 ^, the coefficients t n are 
insensitive to L. However, a very interesting information 
about the development of localization comes from the 
analysis of t n (V) for small samples with L < L c . 

The latter situation can occur in the vicinity of a 
localization-delocalization transition, when the length 
L c (e) can be larger than L. In this critical regime, the 
multifractality is a manifestation of pre-localization in a 
piece of a matter which still dominantly shows the metal- 
lic properties. In the language of the coefficients t n) this 
is described as 



Vt n (e) oc , r (n) = (n - l)d* (n) 



(3) 



where d* (n) differs from the physical dimension d of the 
system and is a function of n. This function gives the 
values of the fractal dimensions d* (n) for each n. For 
ballistic waves in a billiard, one would have d* (n) = d. 
Once we assume that the envelope of a typical wave func- 
tion at the length scale shorter than L c obeys a power 
law V'( r ) r ~ M with a single fixed exponent /i < d/2 
§,§, the set of IPN's shows d* = d- 2[i different from d 
but the same for all n > dj (2/z). This is when one speaks 
Q of a fractal behavior with the fractal dimension d* . 

If d* (n) is not a constant, this signals about some more 
sophisticated structure of wave functions. They can be 
imagined as splashes of multiply interfering waves at dif- 
ferent scales and with various amplitudes, and possibly, 
with a self-similarity characterized by a relation between 
the amplitude of a local splash of a wave function, t, 
and the exponent fx(t) of the envelope of its extended 
power-law tail. The spatial structure of wave functions 
in the vicinity of the localization transition was a subject 
of extensive numerical studies [||,[|,D • To distinguish be- 
tween the fractal and multifractal behavior in numerics, 
one has to detect and compare the values d*(n) coming 
from, at least, several lowest IPN's. In the pioneering 
works [||, the assertion about the fractality at criticality 
has been made on the basis of studies of d* (2), and a 
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non-trivial dependence of dimensionality d* (2) on disor- 
der in 2D metals Q has been observed. Later on, nu- 
merical study of single-electron states at the transition 
between the plateaus in the regime of the Quantum Hall 
Effect (QHE) 0] has established different fractal d*{n) 
for different IPN's, which supported the idea about the 
multifractality in the critical regime. 

The goal of this Letter is to demonstrate that the mul- 
tifractality of wave functions is a generic property of 2D 
disordered systems as soon as the sample size L does 
not exceed the formal localization length L c . In some 
sense, the two-dimensional system is critical not only in 
the QHE regime, and we consider the limit of relatively 
weak magnetic fields. If disorder is weak (the diffusion 
coefficient D is large, so that 2ttvD ^> 1), the localiza- 
tion length L c is exponentially large and, therefore, the 
system is an ideal object to observe the multifractality at 
a large scale. To anticipate a little, our calculation gives 



d*(n) =2-n(4n 2 vD) 



(4) 



which can be associated with a power-law tail of the en- 
velope with the exponent < 1. The result of Eq. (^) 
is essentially non-perturbative and cannot be obtained 
by any expansion in diffusion modes usual for describing 
weak localization effects. 

To get this result, we use below the supersymmetry 
technique 0. Analytical study of wave functions with 
this method has started recently ||. It was discovered 
that the distribution of local amplitudes of waves in a 
very small - zero-dimensional (OD) - disordered conduc- 
tor or a chaotic cavity obeyed the same laws as those 
predicted in the random matrix theories || . One of the 
achievements of the method is that it makes possible to 
calculate the distribution function of local amplitudes 
f(t) as a whole 1 10 , rather than to reconstruct || it 
from the full set of IPN's. This proves to be very helpful 
for calculations beyond OD. 

Using the supersymmetry technique, one reduces the 
calculation of the distribution function / (t) to evaluation 
of a functional integral with the free energy functional 

F[Q] = / Str [D(X7Q (r)) 2 - 2 7 AQ (r)] dr, (5) 

where D is the diffusion coefficient, 7 is a level width,. 
The notations for the supermatrices Q and A, and for 
the supertrace 'Str' are the same as in Ref. Q. A gen- 
eral form of the functional integral valid for an arbitrary 
magnetic field can be found in Ref. |t0|. In what fol- 
lows, we calculate the distribution function / (t) for the 
unitary ensemble, although the results for the orthogonal 
and symplectic ensembles are similar |lT| . The unitary 
ensemble is technically most simple. At the same time, 
it can be a quite interesting one, due to its relation fl^] 
to the recently suggested description of some strongly 
correlated systems, such as high T c materials or 2D elec- 
tron gases at even-denominator filling factors, in terms 
of fluctuating guage fields |TJ,|lJ]. 



For the unitary ensemble, the function / (t) can be 
equivalently rewritten in a somewhat simpler form 



f(t) = ^mj DQ J ^Str(4 X) Q(r)) 

x5 (t- ^2Str(^ 2) Q(r )) exp (-F [Q]) (6) 
where the tt, ' 1 select from Q its boson-boson sector, 
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and to,3 are the Pauli matrices. The limit 7 — > in Eq. 
@ corresponds to a closed system with elastic scattering 
only, and one could, first, calculate the integral over Q(r) 
and take the limit 7 — > at the end. However, it is 
better to get rid of the parameter 7 at an earlier stage 
by integrating, first, over the zero space harmonics of Q. 

If we restrict ourselves with the integration over this 
harmonics only, we arrive at the OD result of Ref. [fl0|| . 
Nevertheless, this would not be enough for probing the 
states which are the precursors of localization, since 
the latter is related to the non-homogeneous fluctua- 
tions Q(r) 0. Moreover, as has been recently shown 
by Muzykantskii and Khmelnitskii jt5|], the treatement 
of rare localization events in the metallic regime (such as 
long-living states in an open mesoscopic conductor) can 
be advanced using the saddle-point solutions of the super- 
symmetric field theory. Thus, we also account for spatial 
variations of Q non-perturbatively, using the saddle-point 
method like in Ref. [jl5l , so that the integration over the 
zero harmonics is only our first step. 

To realize this plan, we represent the matrix Q as 

Q(v) = V(v)AV(v), V(r)V(r) = l, 

and substitute V (r) — » V (r ) V (r) and Q(r) — > 
V (r Q ) Q (r ) V (r D ). In terms of new variables V (r) and 
Q(r), the first (gradient) term in Eq. (||) preserves its 
form, but now the condition Q(r ) = A has to be ful- 
filled. The second term in Eq. (0) transforms into 



F 9 



drStr(Q Q(r)), Q = V (r ) AV (r ) . 



A corresponding substitution can be done in the pre- 
exponential of Eq. (^), too. This enables us to inte- 
grate over Q a . The limit 7 — + drastically simplifies the 
computation, since it makes the essential values of the 
variable 0\ o parametrizing the non-compact sector of the 
supermatrix Q a as large as cosh#i ~ I/7. Performing a 
standard integration over the elements of the superma- 
trix Q - the procedure being a slight modification of the 
calculations of Ref. [|l0| - we reduce Eqs. (||, ||) to the 
form 

/(*) = V& {jQ (ro )=A ex P (-FbQ]) DQ(r)} (7) 
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where the free energy functional F [t, Q] is determined by 
F[t, Q]= J drStr (2|2(VQ (r)) 2 - fAIIQ(r)) . (8) 

The matrix II selects from Q its non-compact boson- 
boson sector Qb = V b AV b : 

TT f K b TT b \ ( WK b 6l/2\ , . 

n= U n)' ViSsatp \« + nBi/2 o J (9) 

and u = cxp(iX7"3), where 6>i > and < x < 27r. The 
OD result can be obtained just by putting Q (r) = A for 
all r, and one ends with the Porter-Thomas distribution 
= l/exp(— Vt) for the unitary ensemble which 
gives the IPN which obey Eq. and show neither mul- 
tifractal nor fractal behavior. 

Nonetheless, this would be only an approximate pro- 
cedure, because the value Q (r) = A does not correspond 
to the minimum of the functional F [t, Q] when t ^ 0. 
The second term acts on the matrix Q as if an external 
field tended to " align" the matrix Q b along such a direc- 
tion that 8i — * oo. However, the condition Q (r G ) = A 
prevents from that, and, as in Ref. |lq| , the minimum 
corresponds to a non-homogeneous configuration of a fi- 
nite Qb (r). 

The solution 9 t corresponding to the minimum of the 
functional F [t, Q] can be found after substituting Eqs. 
(g) into (§) and varying 9\ under the condition nV^f = 
at the boundary of a sample, and 8 t (r Q ) = in the origin. 
As a result, we obtain 

M t(r) = -^j cxp {-6 t (r)) , X (r) = n, (10) 

where A is the Laplacian. This saddle-point equation is 
somewhat different to that discussed in Ref. fl5|] . The 
solution of Eq. (|icj) has to be substituted into the free 
energy functional F [t, Q] which takes the form 

Ft = f{^(\76 t f + texp(-8 t ))dr (11) 

and gives the final result with an exponential accuracy. 

In principle, it can be important to take into account 
fluctuations in the vicinity of the extremum. In the 
lowest order, these fluctuations contribute to the pre- 
exponential. As concerns higher order terms, one should 
estimate them to be sure about the corrections to the ap- 
proximations we made. This can be done in an invariant 
form representing Q (r) as 

W = WA^W.P^ o) (12) 

The supermatrix Vt (r) is determined by Eq. (||) with 
the optimal solution t substituted as the variable 6\, 
and P is treated as a fluctuation around the inhomoge- 
neous optimal solution. To evaluate the contribution of 
fluctuations P, one has to substitute Q (r) from Eq. Jl~2] ) 
into Eqs. ([?], ||) and, then, expand them in P. In the 
zeroth order in P, Eq. (|l2|) corresponds to the optimal 



solution which defines Ft , Eq. ( |il"| ) . After expanding the 
action F[t, Q] in series on the perturbation P, one can 
check that the above choice of Vt cancels the linear terms 
pW, and F[t,Q] = P t +p( 2 ) +p( 3 ) +p( 4 ) + .... At this 
stage of calculating the integral in Eq. (R), one has to 
keep both P t and F^ in the exponent e~ F ^'^ , but ex- 
pand it with respect to other higher order terms on P. 
Then each term of this expansion is given by a Gaussian 
integral and can be evaluated using Wick's theorem. 

Estimating the contribution of the higher-order terms, 
one can see that they are not important as long as the 
sample size L ~ V x l d is smaller than the localization 
length L c which holds in ID and 2D samples. As con- 
cerns the 3D case, the fluctuations are always small when 
disorder is weak. In all the cases, the contribution from 
F t dominates, but we are also aware of what comes from 
square terms F^ 2 \ Calculating Gaussian integrals one 
obtains for the pre-exponential J 

J = cxp {i £„ In (xj(n)/ U a ,p Xi(«, A «)) } • (13) 

J differs from unity, since the optimal fluctuation breaks 
the symmetry between anti- (f) and commuting (b) 
components of P and splits the spectra {Xf( n )} an d 
{Xb(oi, /3, n)} of their eigen-modes. The set of four- 
fold degenerate modes Xf an d non-degenerate Xb are 
defined from the 'Schroedinger equations' [— 2-kvD/S. + 
U — x] s = 0, where the effective potentials P(/.b) Q = 
2itvDk( f.bkfV^) + t K (f.b),a ex P( — 8t) {k and k are real 
numbers pi] ) obey the sum rule U b ,a = 4t//>. 

Until now all the manipulations did not depend on the 
dimensionality. However, solving Eq. ( ^0| ) has to be done 
for each dimensionality separately. Results for ID and 
3D will be presented elsewhere pT| . Below, we discuss 
the 2D case. The use of the a-model limits the length- 
scales from below by the value of the mean free path /. 
Therefore, when solving Eq. (|i~0|), we cut off the radii 
r < I and replace the conditions Q = A (dt = 0) at the 
origin by the same requirement at r — ro ~ I, provided 
the dependence of the finally derived optimal action F t 
on the cut-off length is weak. For the sake of simplicity, 
we consider the sample in the form of a disk of a radius 
L, ln(L/Z) (2tti/D) 2 , and place the observation point 
r D in its center. This geometry allows us to seek for an 
axially symmetric solution of Eq. ( |i~0| ) which, with the 
parameter p = \J2ttvD jtl 2 , takes the exact form 

where A — 1 — fi is determined from the boundary con- 
ditions as 

V^Tp^+A = ^f^. (14) 

Eq. (|l4|) has positive roots if p > In j 3> 1, which limits 
the amplitude t by the value (XfI)^ 1 , the density of a 
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plane wave forward-and-backward scattered within the 
mean free path length. In this limit ji -C 1, which is 
the reason for the power-law asymptotics of the optimal 
exp(— 6t). The same conditions also simplify the search 
of roots of Eq. ( |l4| ) and reduce it to the equation on fr. 



z(T) 
V ~ 21n(L/i) ' 



? p z _ T = tVln(L/l) 
ZL — 1 — 2n->uD ■ 



(15) 



The leaading terms of the optimal free energy F t can 
be found as 

F t «(2 7 r) 2 ^{ A1 + M 2 lny}, 

and the logarithmic dependence of F t on the ratio L/l 
in all the cases justifies the cut-off procedure and makes 
meaningful the use of the derived formulae for an arbi- 
trary position of the observation point and the shape of 
the sample. 

One can also evaluate the pre-factor J from Eq. (|l£ 
and find that J = 1 + 0(T 2 ) at T < 1, and J oc 
/i exp (/i In j) ~ TlnT when T » 1. This confirms our 
previous statement that the optimal action dominates 
over other contributions and enables us to find the form 
of the distribution function f^ 2 \t) for small, t = \ip\ 2 <, 
and for large, t > 2ir 2 vD/(V\nj), amplitudes with the 
exponential accuracy: 



exp(-Vt[l-T/2]) ; 



exp | 



' \n(L/l) 



ln 2 T 



T< 1 
T> 1. 



(16) 



From this expression, one can conclude that disorder 
makes the appearence of high-amplitude splashes of wave 
functions much more probable than one would expect 
from the OD Porter-Thomas formula. Obviously, this 
is a localization effect. But the tails of the states do 
not decay exponentially, as one would expect for a par- 
ticle localized in a quantum well or in a ID wire beyond 
the localization length: Even in the asymptotic regime, 
the size L of the system influences the distribution. The 
splashes look as if they were formed by focusing the waves 
by those rare configurations of scatterers which play the 
role of 'parabolic' mirrors. The structure of these states 
can be anticipated from the way how their distribution 
function feel the boundary. Since this comes in through 
the form of the optimal fluctuation, we believe that they 



have the power-decaying tails, |V>( r )l^ 



(l/r) 



2p 



(/i < 1), which approach the limiting r~ 2 dependence for 
the highest amplitudes t ~ (IXf)^ 1 ■ 

Moreover, the form of IPN's, t n derived on the basis of 
Eq. (|ilf ) shows a scaling with the size of a system which 
forces us to assign them a multifractal nature. To find the 
moments t n accurately enough, we have to take into ac- 
count that, despite of that the cross-over to the Od case 
looks like a formal limit T(t) — > 0, the Porter-Thomas 
statistics fails unless the condition tV <C \J2txvD is satis- 
fied. Hence, only the first few ratios t n , 2 < n <C V2ttuD, 
can be estimated by using a finite polynomial expansion 



of f(t) into the series on T. The result would look as 
a correction to the values expected for chaotic quantum 
billiards Alternatively, we derive the higher order 

IPN's from the intermediate result of Eq. ((?]), 



n %n^ J dtJ ^) ex P + (" - 2 ) HtV)} . 



by the saddle-point method. The moments t„ calculated 
in both ways are in a good agreement with each other 
and, in the leading order, take the form 



min{n!, [2nvD/]n k ] n } ( 1 
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V 



i-6 



n(n — 1) 
8n 2 vD ' 



To summarize the result of the derived above ex- 
act statistics of amplitudes of wave functions in meso- 
scopic disordered conductors, we end with such a volume- 
dependence of the inverse participation numbers t n which 
manifests the multifractal behavior of quantum states, 
Eq. (|J). Multifractality seems to be the generic property 
of 2D disordered systems. In the unitary ensemble, the 
multifractal dimensions are d* — 2 — n(47r 2 ^D) _1 , and 
all over the metallic regime, the dependence of d* on n 
and the level of disorder, Eq. (^), is exact and can be a 
subject of numerical verifications. 

We thank D. Khmelnitskii for sending us the preprint 
MM prior to publication. 
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